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NOTES ON GALERKIN-FINITE ELEMENT METHODS FOR THE SHALLOW 
WATER EQUATIONS WITH CHARACTERISTIC BOUNDARY CONDITIONS. 


D.C. ANTONOPOULOS AND V.A. DOUGALIS 


Abstract. We consider the shallow water equations in the supercritical and subcritical cases in one space 
variable, posed in a finite spatial interval with characteristic boundary conditions at the endpoints, which, as 
is well known, are transparent, i.e. allow outgoing waves to exit without generating spurious refiected waves. 
Assuming that the resulting initial-boundary-value problems have smooth solutions, we approximate them 
in space using standard Galerkin-finite element methods and prove error estimates for the semidiscrete 
problems on quasiuniform meshes. We discretize the problems in the temporal variable using an explicit, 
fourth-order accurate Runge-Kutta scheme and check, by means of numerical experiment, that the resulting 
fully discrete schemes have excellent absorption properties. 


1. Introduction 

In this paper we consider the system of shallow water equations 

r]t + Ux + {'nu)x = 0 , 

Ut+Vx + UUx = 0 , 

a well known approximation of the two-dimensional Euler equations of water-wave theory, modelling two-way 
propagation of long surface waves of finite amplitude in a uniform horizontal channel of hnite depth, |Wh) . 
The variables in (dU are non-dimensional and unsealed; a; G R and t > 0 are proportional to position along 
the channel and time, respectively, while rj = r]{x, t) and u = u{x, t) are proportional to the elevation of 
the free surface above a level of rest and to the horizontal velocity of the fluid, respectively. The latter is 
depth-independent to the order of approximation represented by the scaled analog of (ini). In the variables 
of (|l.ll) the bottom of the channel lies at a depth equal to —1. 

It is well known that the initial-value problem for (dH), posed with smooth initial data r]{x, 0) and m(x, 0) 
for a; S K, has, in general, smooth solutions only locally in t, cf. e.g. [Maj . Ch. 2. In this paper we will pose 
(dH) in the finite ‘channel’ [0, L] with given initial values at t = 0, 

r]{x,0) = r]°{x), u(x,0) = u^(x), 0<x<L, (1.2) 

and consider transparent boundary conditions at x = 0 and x = L, i.e. conditions that permit the waves to 
exit the ‘computational’ domain [0, L] without generating spurious reflected waves that pollute the solution 
inside [0,L]. The transparent boundary conditions that we will use are nonlinear characteristic boundary 
conditions for subcritical and supercritical flows governed by (dH. Such conditions were first used, to our 
knowledge, by Nycander et al., [NMFj . in numerical experiments with finite difference discretizations of 
shallow water models. In the paper at hand we will analyze Galerkin-finite element approximations for 
smooth solutions of the initial-boundary-value problems (ibvp’s) resulting from the application of charac¬ 
teristic boundary conditions to mi). The well-posedness of these ibvp’s was analyzed in Petcu & Temam, 
[PTl] and Huang et al., [HPTj . 

The characteristic boundary conditions may be derived as follows. We write the system mu as 
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^ y / eigenvalues Ai = u + -^1 + 77, A2 = m — -y/l + 77. Assuming always that 
77 > —1, we consider two types of flows: 

Supercritical : u > + 77, 0 < A2 < Ai, 

Subcritical : u < \/l + 77, A2 < 0 < Ai. 


It is well known , cf. e.g. |Wh| . that along the family of characteristic curves with x(t) = Ai = u-\-yJ\ + 77, the 
quantity ri := m + 2i/1 + 77 is constant, while along the curves with x{t) = A2 = u — ^/T^^rj, r2 ■= u — 2^/1 + 77 
is preserved. If 77 and u are expressed in terms of the Riemann invariants ri, r2 one obtains the system of 
ordinary differential equations (ode’s) 


dri „ , , 

-j- = 0 on x{t) : 

= 0 on x{t) : 


^i=Mrur.) 

|=A2(r„r2) 


3ri + r2 
4 

ri + 3r2 
4 


which is equivalent to the original pde system (ED and whose discretization yields the classical method of 
characteristics for solving (HID- If we pose (ED in the spatial interval [0,L] it is straightforword to see, 
cf. [Whj . Section 5 . 4 , that the temporal integration of the ode system requires in the supercritical case that 
ri( 0 ,t) and r2(0,t) be given for t > 0 , as both families of characteristics are incoming at a; = 0 ; this is 
equivalent to prescribing u{ 0 ,t) and 77(0, t) for t > 0 . In the subcritical case 7 'i( 0 ,t) and r2{L,t) should be 
given for t > 0 , since they correspond to the incoming characteristics at a: = 0 and at x = L respectively. 

Following |NMF) we assume that outside the interval [ 0 , L] the flow is uniform and is given by ri{x, t) =770, 
u{x,t) = uq, where 770, uq are known constants. Therefore, in the supercritical case the characteristic 
boundary conditions are simply 

77(0, t) = 770, u{ 0 ,t) = uo, ( 1 . 3 ) 

with Uq > \/l + 770 , while, in the subcritical case, they are of the form 


77(0, t) + 2 -\/l + 77(0, t) — uo + 2 -\/l + 770, 

_ _ ( 1 . 4 ) 

u{L, t) - 2 ^J\ + ri{L,t) = uo- + 770, 

where now it is assumed that itQ < 1 + 779. In both cases we may view the solution (77,77) of (jl.ll) for 
0 < X < L, t > 0 , generated by the initial conditions ED, as a perturbation of the uniform flow (770,779) 
to which the solution inside the computational domain [0,L] will revert once the waves generated by the 
initial conditions exit this interval. It is straightforward to check, using the definitions of characteristics and 
Riemann invariants and considering e.g. initial conditions that differ from 779, uq in a subinterval of [ 0 ,L], 
that the boundary conditions ED and (ED are transparent. 

As previously mentioned, the characteristic boundary conditions ED and (ED were used by Nycander 
et ai, |NMFj . in finite difference simulations of the shallow water equations, in one space dimension, actually 
in more complicated instances of hydraulic and geophysical interest, including single- and two-layer flows 
in channels of variable width and variable bottom topography, time-dependent forcing in the boundary 
conditions, examples where transcritical flows develop, et al. (In the case of two-layer flows an approximate 
SW system was used in which the barotropic and baroclinic modes are decoupled; this allows using the 
analogs of the (local) characteristic boundary conditions in this case too. Also, for reasons of numerical 
stability, a diffusive term with a small viscosity coefficient was added in the momentum equations.) The 
finite difference spatial discretization was effected on a staggered grid and the leap-frog scheme was used for 
time stepping. Similar model equations and characteristic boundary conditions were applied in simulations 
of two-layer hydraulic exchange flows in |FM) . 

The characteristic boundary conditions ED and (ED were also used by Shiue et al. , |SLTTj , in the case 
of the one-dimensional, single-layer shallow water equations in channels of variable bottom topography in 
the presence of Coriolis terms and with the addition of a cross-velocity variable that depends only on x. The 
system, written in balance law form, was discretized in space using midpoint quadrature for the source cell 
integral and a ‘central-upwind’, |KNP) . [KPj . Godunov-type approximation of the flux term; a second-order, 
explicit Runge-Kutta method was used for time stepping. Many numerical experiments performed with this 
scheme are reported in [SLTT] : they simulate interesting cases of subcritical, transcritical and supercritical 
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flows over variable-bottom topographies. The characteristic boundary conditions and the same numerical 
scheme were subsequently used in [BPSTT] in the case of two-layer problems in one dimension under the 
decoupling assumptions of |NMFj . (The local well-posedness of this two-layer problem was studied in |PT2) .i 

In case the elevation of the free surface ry is a small perturbation of the steady state rjQ one may derive 
linearized approximations to the characteristic boundary conditions (HI) in the subcritical case. These 
linearized conditions are also considered in |NMF] and in |SLTT] , where they are compared to the nonlinear 
exact conditions and found in general to cause spurious reflections that enter the computational domain. 
(The linearized boundary conditions are easily seem to be (exactly) transparent for the linearized shallow 
water equations obtained by linearizing (ED about the steady state {r]o,uo). In [SLTTj it is shown that the 
ibvp for the linearized system supplemented by the linearized boundary conditions is well posed.) As pointed 
out in [NMF] . [SLTTj . and in m, the linearized characteristic boundary conditions have been extensively 
used in the computational fluid dynamics literature; the last reference contains a review of several other 
absorbing boundary conditions for the shallow water equations at artificial boundaries, including absorbing 
(‘sponge’) layer conditions et ai. 

Of particular interest for our purposes is the rigorous analysis of the well-posedness of the ibvp’s (locally 
in time) for the shallow water equations with characteristic boundary conditions carried out in Huang et 
ai, |HPT] . and Petcu & Temam, [PTTj . The ibvp in the supercritical case, was studied by Huang et ai, 
[HPT] , in fact in the more general setting of shallow water supercritical flows over a variable bottom in the 
presence of Coriolis terms and a lateral component of the horizontal velocity depending on x, and also with 
nonhomogeneous boundary conditions satisfying appropriate compatibility conditions. The hypotheses of 
[HPT] on uq, rjo and the initial data, briefly reviewed in section 2 in the sequel, guarantee the existence 
and uniqueness, locally in time, of a smooth solution of the ibvp (ll.ll) - (ll.3l) with positive 1 + rj, satisfying 
the strong supercriticality property — (1 -I- 77 ) > Cg for some positive constant cg. The well-posedness 
of the ibvp in the subcritical case, i.e. of the ibvp (ED, (ED, (ED, was studied by Petcu & Temam, 
[PTlj . The assumptions of [PTlj (reviewed in section 3 below) imply the existence and uniqueness, locally 
in time, of a smooth solution of the ibvp with positive 1 -I- 77 and satisfying the strong subcriticality condition 
77 ,^ — (1 -I- 77 ) < —Cg, where cg is a positive constant. 

In this paper we will analyze standard Galerkin-finite element spatial discretizations of the ibvp ’s (ED- 
(ED and (fTTI) . (fm . (fLl) . under the hypothesis that they have smooth solutions. In both cases the basic 
approximation will be effected by functions which are piecewise polynomials of degree r — 1, r > 2 , 

on quasiuniform partitions of [0,L]. In section 2 we consider the supercritical case and prove an L^-error 
estimate of accuracy for the Galerkin approximations of 77 and u. (It is well known that this is the 

expected best order of convergence in for standard Galerkin semidiscretizations of first-order hyperbolic 
problems on general quasiuniform meshes. For uniform meshes, better results hold, cf. m for the analysis 
in the case of a linear model problem. In |AD2j it was proved that the order of convergence in for 
piecewise linear continuous elements on a uniform mesh is equal to 2 in the case of an ibvp for (ED with 
the homogeneous boundary conditions u{0,t) = u{L,t) = 0. This superaccuracy result is expected to hold 
for the ibvp’s under consideration as well and this is indeed what the numerical experiments of section 4 
indicate.) For the proof of the error estimate we assume that a strengthened supercriticality condition holds 
for the solution of (irT])-(f 01 ): cf. (H1)-(H3) in section 2. The proof also requires that r > 3 so that a certain 
bootstrap argument, based on the boundedness of the || • ||i_oo norm of an error term, goes through as in 

[mj . [AD^ . 

In section 3 we turn to the subcritical case. We write the ibvp (fTTI) . dEH), (HI) in its classical diagonal 
form in which the new unknowns are analogs of the two Riemann invariants in the context of the ibvp 
at hand and satisfy homogeneous Dirichlet boundary conditions one at x = 0 and the other at x = L. 
The diagonal system is discretized in space on a quasiuniform mesh by the same type of standard Galerkin 
method as before, and a error estimate of is proved for both components of the solution. A 

change of variables of this semidiscrete approximation yields approximations of the original unknowns 77 and 
u of accuracy. The proof requires that a strengthened form of the subcriticality property holds for 

the solution of the ibvp (11.11) . (II.2|) . (11.41) (cf. (lYll) . (IY2|) in section 3), and the technical assumption that 
r > 3. 

Section 4 is a report of various numerical experiments that we performed with the Galerkin-finite element 
methods of sections 2 and 3 and some of their variants. We use spatial discretizations with piecewise linear 


3 




































continuous functions on uniform meshes and discretize them in the temporal variable by the ‘classical’, 
explicit, four-stage, fourth-order Runge-Kutta scheme. The resulting fully discrete methods are stable under 
a Courant number restriction. (Stability and convergence of high order explicit Runge-Kutta methods 
was established for closely related pde systems in |AD1] and |AD2] .l Our main purpose in the numerical 
experiments is to check the stability and the numerical order of convergence of the fully discrete Galerkin 
methods and study by computational means their absorption properties. Although the full discretizations 
of the characteristic boundary conditions are not exactly transparent of course, the numerical experiments 
show that they are practically transparent, in contrast to the analogous Galerkin schemes with the linearized 
boundary conditions that we also implement in the subcritical case; the latter are absorbing but in general 
allow spurious reflections to form and enter the computational domain. In the subcritical case we also 
implement the analogous fully discrete Galerkin method for the original, nondiagonal form (dl), (HH), 
dUD of the system and check that it gives results close but somewhat inferior to those of the analogous 
discretization of the diagonal form of the system. 

In the error estimates in the sequel, we let the spatial interval be [0,1] for simplicity. We let = C^[0, 1], 
fc = 0,1, 2,..., be the space of k times continuously differentiable functions on [0,1]. The norm and inner 
product on = L^(0,1) are denoted by || • ||, (•,•), respectively. For integer fc > 0, || • ||fc will 

denote the usual, L^-based Sobolev spaces of classes of functions and the associated norms. The norms on 
L°° = L°°(0,1) and on the L‘^-based Sobolev spaces = W^{0, 1) will be denoted by || • ||oo, || ■ |U,oo, 
respectively. Finally, we let be the polynomials of degree < r. 


2. Semidiscretization of the supercritical shallow water equations 


In this section we consider the shallow water equations with characteristic boundary conditions in the 
supercritical case. Specifically, for {x, t) G [0,1] x [0, T] we seek p = r]{x, t) and u = u(x, t) satisfying the 
ibvp 


+ Mx + (lyR)x = 0, 

Wt + + UUx = 0, 

r]{x,Q) = rf{x)^ u(a;, 0) = 
r]{0,t) = T]Q, u{0,t) = uo, 


0<a;<l, 0<t<T, 

= u°{x), 0 < a: < I, 

0<t<T, 


(SWI) 


where 7?°, are given functions on [0,1] and tjq, uq constants such that I + > 0, Uq > 0, Uq > \/l + pg. 

As mentioned in the Introduction, the ibvp (ISWIj) was studied by Huang et al., [HPT], in fact in the 
more general case of a shallow water supercritical flow with nonhomogeneous boundary conditions over a 
variable bottom for a nonzero Goriolis parameter and also in the presence of a lateral component of the 
horizontal velocity depending on x only. In the simpler case of (ISWII) , the proof of the main result of |HPT| 
amounts to the selection of a suitable constant solution ( 770 , wq) of (ISWII) and of sufficiently smooth initial 
conditions close to the constant solution and satisfying appropriate compatibility relations at a: = 0. Under 
these hypotheses the conclusion of [HPT] is that given positive constants cq, ag, and Coi there exists a 
T > 0 such that a sufficiently smooth solution of (ISWII) exists satisfying for {x,t) G [0,1] x [0,T] the strong 
supercriticality properties 

U^-{1 + V)>cl (PI) 

u > ag, (P2) 

Co<(l + ^)<Co- (P3) 

For the purposes of the error estimation to follow we will assume that (ISWII) has a sufficiently smooth 
solution (ry, u) that satisfies a strengthened supercriticality condition of the following form: There exist 
positive constants a and /3, such that for {x,t) G [0,1] x [0,T] it holds that 


1 +77 > /3, 
u > 2a, 

1 + 77 < {u — a){u — ^). 
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(HI) 

(H2) 

(H3) 























Obviously m and (IH3I) imply that u > y/1 + rj. It is not hard to see that (IH3I) follows from (IP1I) - (IP3I) if 
e.g. Qfo is taken sufficiently small and cq sufficiently large. We also remark here that in the error estimates 
to follow (IH3I) will be needed only at cc = 1 for t G [0, T], 

We will approximate the solution of (ISWII) in a slightly tranformed form. We let rj = rj — tjq, u = u — uq 
and rewrite (ISWip as an ibvp for rj and u with homogeneous boundary conditions. Dropping the tildes we 
obtain the system 


It + Moba; + (1 + VojUx + {r]u)x = 0, 

„ 0 < a: < 

Ut + r]x + uqUx + uUx = 0, 

ri{x, 0) = rf{x) - rio, u{x, 0) = u°{x) - uo, 


1, 0 < t < T, 
0 < a; < 1, 


(SWla) 


77 ( 0 , t) = 0, u{0, t) = 0, 0<t<T. 


In terms of the new variables 77 and u, our hypotheses (lHT])-(lH3l) become 


l + V + Vo^ 13, (Hla) 

M + 77o > 2a, (H2a) 

1 + 77 + 770 < (m + Mo — ce){u + Mo — ^). (H3a) 


In the sequel, for integer fc > 0, let (7^ = {m G C^[Q, 1] : 7;(0) = 0}, and = {^ G 1) : v{0) = 

0}. For a positive integer let 0 = xi < a ;2 < • • • < Xat+i = 1 be a quasiuniform partition of [0,1] with 
h := maxAxij^i — xA, and for integer r > 2 define Sh = {(/> & , G Pr-i, 1 < j < .^V}. It is 

I [Xj ,Xj + i \ 

o o 

well known that if m G , there exists x & Sh such that 


\\v-x\\+h\\v'-x'\\<Ch^\\v\\r, (2.1) 

and, cf. |Sch] . if r > 3, 

\\v-xh<Ch^-^\\v\\r. (2.2) 

(Here and in the sequel C will denote a generic constant independent of h.) In addition, if P is the L^- 
projection operator onto Sh, then it follows that, cf. [dTTw] . 


IIFmIIoo < CIImIIoo, if mGL“, (2.3) 

\\Pv-v\\^<Ch^\\v\\r,^, if (2.4) 

As a consequence of the quasiuniformity of the mesh the inverse inequalities 

llxlli<C/7-i||x||, (2.5) 

||xIU<C/7-(^+i/2)||;^||, _^- = o,1, (2.6) 

hold for X G Sh- (In dM]) || • ||o,oo = || • lU-) 

The standard Galerkin semidiscretization of (ISWlal) is defined as follows: We seek r]h, Uh ■ [0,T] —>• Sh 
such that for 0 < t < r 

o 

{'nht,(l>) + {uo'nhx,4') + {{i + m)uhx,(t)) + {{'nhUh)x,4') = o, yr/rGSh, ( 2 . 7 ) 

o 

{Uht,<l)) + {Vhx,<l)) + {uoUhx,<l}) + {uhUhx,(t}) = 0, y(t)&Sh, ( 2 . 8 ) 


with 

?77i(0) = F(77°(-) - 770 ), Mft,(0) = P(m°(-) - Mo). (2.9) 

The main result of this section is: 


Proposition 2.1. Let (rj,u) he the solution of iSWld) . and assume that the hypotheses \Hla\) . !lH2a\} . /filial) 
hold, that r > 3, and h is sufficiently small. Then the semidiscrete ivp |^. 7| )- (KB has a unique solution 
{r]h,Uh) for 0 <t <T satisfying 

max (||77(t) - r]h{t)\\ + \\u{t) - Uh{t)\\) < Ch^~^. (2.10) 
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Proof. Let p = rj — Pp, 9 = Prj — a = u — Pu, f, = Pu — Uh- After choosing a basis for Sh, it is 
straightforward to see that the semidiscrete problem represents an ivp for an ode system which 

has a unique solution locally in time. While this solution exists, it follows from (lO) . (12:81) and the pde’s in 
(ISWla|) . that 

{6t, (ft) + {uo{px + 9x), 0) + ((1 + »7o)(o'i + fx), (f) + {{r)u - T]hUh)x, (f>) = 0, V<?!) £ Sh, (2-11) 

[ftA) + {px + dx,(j)) + {uo{(Tx + £,x),(l>) + {uux - UhUhxA) = y(t>€Sh. ( 2 . 12 ) 

Since rju - phUh = p[a- + ^) + u{p + 6») - (p + e){a + ^), uu^ - uhUhx = {u(t)x + - crcTx - ^^x, 

it follows that 


{pu - phUh)x = {pOx + {u0)x - i9^)x + Ri, 

UUx - UhUhx = {uOx - ^£.x + R 2 , 

where 

Ri = {w)x + {up)x - {p(^)x - {pCjx - {Q(^)x, 

R 2 = {ua)x - {(J^)x - crdx- 

Therefore, the equations (12.111) . (12.121) may be written as 

{Ot, 0) + {uoex,f>) + {'fxA) + {{ue)x, <t>) - {{eOx, <t>) = -{RiA), e Sh, 

(6, + {dx, <!>) + {uofx, 0) + («).. - i^fx, 0) = -{R2,f>), V0 € Sh. 

where 7 = (1 + 770 + p)f. and 

Ri = uqPx + (1 + Po)(Jx + Ri, 

R2 = Px + Uoax + R2- 

Putting if = 6 in (I2.17p . using integration by parts, and suppressing the dependence on t we have 
kimf - (7, 9 x) + i(uo + u(l))e^(l) + (1 + 70 + 7 ( 1 ))C( 1 ) 0 ( 1 ) 

- ic(i)e'(i) = -^{Ux9, e) + e) - (i?i, e). 

Take now f = Pi = ^[(1 + Vo + v)^] in (12.181) and get 

(6,7) + { 9 x,l) + {uofx,l) + («)^, 7) - (CCx, 7) = -(A3, P7 - 7) - (A2, Pp), 

where 

R 3 — ^x P '^o^x “f {'^^)x red¬ 
integration by parts in various terms in (I2.22p gives 

{uo^x,l) = {uo^x, (1 + 7 ?o + v)0 = 5 ^ 0(1 + 7 o + »?(l))r^(l) - ^{uoVx^, Oi 

(«)d,7) = («)d, (1 + 7?o + V)0 = {ux^i (1 + % + V)0 + «d, (1 + 70 + V)0 

= 5 “(i)(i + 70 + 7 (i))r^(i) + i(Ud(i + 70 + 7 ), r^) - U^vxf,, r), 

(rrd,7) = (rrd,(i + 7o + 7)r) = i(i +70 + 7(i))r^(i) - i(7dr^r)- 

Hence (|2.22ll becomes 

( 6 , 7 ) + {Sx,l) + i(Mo +w(l))(l + 70 +7(l))r^(l) - ^(1 + 70 + 7(l))r^(l) 

= (A4,r)-(A3,P7-7)-(A2,P7), 

where 

R 4 = ^UoPx^ - ^Uxil + 7o + 7)r + \upxi - \vx^'^- 
Adding now (12.211) and (12.241) we obtain 

^Jj[ll^lP + ((l + 7o + 7)^1 r)] + = ^{pt^^O — S) 

+ U^xO, 9) - (Pi, 7) + (P 4 , r) - (P 3 , P 7 - 7 ) - {R 2 .P 1 ). 
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(2.13) 

(2.14) 

(2.15) 

(2.16) 

(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 


(2.24) 

(2.25) 


(2.26) 




















where 


UJ =i(Mo + u(l))0^(l) + i(Mo + m(1))(1 +770 + 77(1))C^(1) 

+ (1 + 770 + 77(1))C(1)0(1) - - 5(1 + ^0 + v(l))ea)- 

In view of (12.9L by continuity we conclude that there exists a maximal temporal instance > 0 such that 
iVh,Uh) exist and ||Ca:||oo < a for t < th- Suppose that th < T. Then, since ||^||oo < ||^a;||oo, it follows from 
(12.271) that for t G [0, th] 


w >5(770 + 77 ( 1 ) - a ) 9 ‘^{ l ) + i(l + 770 + 77 ( 1 ))( 77 o + 77(1) - 
+ (l + 77o + 77(lMl)0(l) = K0(l),m)f 


(2.28) 


where fx = uq + 77(1) — a, A = 1 + 770 + 77(1), ix = uq + 77(1) — The hypotheses (IHla|) and (IH2a|) give that 

0 < /7 < 77, A > 0. It is easy to see then that the matrix in (12.281) will be positive semidefinite precisely when 
(IH3al) holds. We conclude from (I2.28|) that w > 0. 

We now estimate the various terms in the right-hand side of (12.261) for 0 < t < th- We obviously have 

l(r77?,C)l + l(«.0,0)l<C(||ef+ ||0f), (2.29) 


1(^x0,0)1 <a||0f. (2.30) 

In addition, from (|2.19|) , (I2.14|) , and the inverse and approximation properties of Sh and (12.31) , (12.41) we have 

\iRi,o)\ <ch^-^\m + iipiiooiiexiiii0ii + iipx||oo||?iiii0|| 

+ lkx||oo||0f +|k||oo||0x|l||0|| (2.31) 

<c/7’-i||0||+c(ii0f+ iief). 


Also, from (|2.25l) 

\{R 4 ,o\<cm^+cmuur<c{i+a)m^. 

By (11201), dug), (ESI), (EH) and the inverse and approximation properties of Sh we have 


\{R 2 ,Pi)\<ch^-^\m+c\\a\u\a\m+ch^^^ 

<ch^-^m+cur. 


(2.32) 


(2.33) 


Finally, using a well-known superapproximation property of Sh^ cf. [rTnw] . [D 2 ], in order to estimate the 
term P 7 — 7 by 

11^7 - 7ll = ll^[(l + ^ + 7o)C] - (1 + 7 + Vo)a < ChUW, (2.34) 

we obtain by (12.231) and the inverse properties of Sh that 

|(i?3,P7-7)| < |(6I,,,P7-7)| + \{uo^x, P'1 - 7)\ + l(«)x,-P 7 - 7)1 + l(??x,-P 7 “ 7)1 

< Ch\\ex\m\\ + C/7||?x||||?|| + ChUr + C/7||?||oo||?x||||CII 

^C'll^llllCII+qiCf+ Ca||Cf ^ 

<C||0f+ C(l + a)||ef. 

Therefore, (|2.26l) . the fact that w > 0, and the inequalities (I2.30p - (I2.33L (12.351) give ior 0 < t < th 

iiPf + ((1 + ^0 + 7)e,C)] < Ch^-\m + lieil) + C{\\9f + IlCf), 

where C is a constant independent of h and th- By (IHlal) the norm ((1 + 770 + 77 )-, •)^/^ is equivalent to that 
of uniformly for t G [OiT]. Hence, Gronwall’s inequality and the fact that 0(0) = ^(0) = 0 yield for a 
constant C = C{T) 

||0|| + ll^ll < for 0<t<th- (2.36) 

We conclude from (12.61) that ||Cx||oo < ChT~^l'^ for 0 < t < th, and, since r > 3, if /7 is taken sufficiently 
small, th is not maximal. Hence we may take th = T and (|2.10l) follows from (12.361) . □ 
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3 . Semidiscretization of the subcritical shallow water equations 

In this section we consider the shallow water equations with characteristic boundary conditions in the 
subcritical case. Specifically, for {x,t) € [0,1] x [0,T] we seek rj = ri{x,t) and u = u{x,t) satisfying the ibvp 

Vt+Ux + ir]u)x = 0 , 

Ut+r]j:+ UUx = 0 , 

r]{x,Q) = rf{x), u{x,Q) = u^{x), 0 < a; < 1, (SW2) 


0 <a;<l, 0<t<T, 


u(0, t) + 2-\/l + ?7(0, t) — Uq + 21 + ?7o, 0 < t < T, 

u{l,t) - 2^1 + ri{l,t) = uo - + r]o, 0<t<T, 

where 7?°, are given functions on [0,1] and tjq, uq constants such that 1 + t^q > 0 and Uq < 1 + yo- 

As mentioned in the Introduction, the ibvp (jSW2|l was studied by Petcu & Temam, [PTI] . They used the 
hypotheses that there exists a constant co > 0 such that Uq — (I + rjo) < —Cg and that the initial conditions 

7)^(x) and u^(x) are sufficiently smooth and satisfy the condition (u°(x))^ — (l + r/^(x)) < — c§ (with l + 7]^(x) 

positive) and suitable compatibility relations at a; = 0 and a; = 1. Under these assumptions one may infer 
from the theory of [PTI] that there exists a T > 0 such that a sufficiently smooth solution (r/, u) of (ISW2I) 
exists for (x, t) G [0,1] x [0, T] with the properties that 1 + 77 is positive and the strong subcriticality condition 


- (1 + 17 ) < -Co, 


(n) 


holds for (x, t) G [0,1] x [0, Tj. For the purposes of the error estimation to follow we will assume that (ISW2I) 
has a sufficiently smooth solution ( 77 , u) such that 1 + 77 > 0 and satisfies a stronger subcriticality condition. 
Specifically we assume that for some constant cq > 0 it holds that 

Wo + \/l + i?o > Co, Wo - \/l + 770 < -Co, (Yl) 

and for {x,t) G [0,1] x [0,T] that 

u + \/l + 77 > Co, u — -^/l + 77 < —Co. (Y2) 

Obviously (jYIl) and (jY2l) imply the subcriticality conditions Mq — (I + 770 ) < —Cq and (jHI of jPTIj . and 
approximate the latter better as co decreases. 

In this section we will approximate the solution of (ISW2|) with a Galerkin-finite element method after 
transforming (ED in its classical diagonal form. As in the Introduction, we write the system as 


+ A 


= 0 , 


( 3 . 1 ) 


where A = 


The matrix A has the eigenvalues Ai = u + y/1 + 77 , X 2 = u — -v/l + 77 , (note 


' u 1 + 77 ^ 

1 u 

that by (lY2l) Ai > Co and A 2 < —cq in [0,1] x [0,T]), with associated eigenvectors Xi = (-v/l + 77 ,1)^, 
X 2 = (—-^1 + 77 , 1 )^. If S is the matrix with columns Xi, X 2 it follows from (13.1|) that 


S 


-1 (rit 

Ut 


+ 


Ai 

0 


0 

A 2 


S 


-1 1-. 


= 0 . 


(3.2) 


If we try to define now functions u, w on [0,1] x [0,T] by the equations S 


-1 (m 
Ut 


.-1 / Vx 


we see that these equations are consistent and their solutions are given by u = h{u + 2-^1 + 77 ) + c„. 


w = i(u — 2y/l + 77 ) + Cyj, for arbitrary constants c„, c^,. Choosing the constants c„, Cw so that v{0,t) = 0, 
777 ( 1 , t) =0, and using the boundary conditions in (ISW2I1 we get 


77 = i[7t - ito + 2(\/l + 77 - do)], 777 = i[u - 77o - 2 (y^1 + 77 - do)], 

where dp = \/l + rjQ. The original variables 77 , u are given in terms of v and w by the formulas 


(3.3) 


V = [^(w — w) + do]^ — 1, 77 = 77 + 777 + 770. 


(3.4) 



































Since 


Ai — u + \/l + 77 — uo + ^0 + 


3v + w 


\2 = u- v'l + 77 = uo - (5o + 


V + Stc 


(3.5) 


we see that the ibvp (ISW2I) becomes 


0 


0 + ^ 

0 

n\ — „„0 


= 0, 0 < X < 1, 0 < t < T, 


7;(x, 0) = 7;“(x), 7i;(x, 0) = 7c^(x), 0 < x < 1, 

77 ( 0 , t) = 0, 7/;(l, t) = 0, 0 < t < T, 


(SW2a) 


where x°(x) = ^[u°(x) — mq + 2(y^ 1 + 77°(x) — Jq)], 7i’°(x) = \[vP{x) — uq — 2(-\/l + 770 (x) — (5o)]. Under our 
hypotheses (|SW2aP has a unique solution (u, ic) on [0,1] x [0, T] which will be assumed to be smooth enough 
for the purposes of the error estimation that follows. Of course, v and w represent analogs of the Riemann 
invariants of the shallow water system in the context of the ibvp at hand; the system of pde’s in (ISW2all 
and (EH), (EH imply that the solution ( 77 , u) of (ISW2I) may be expressed in terms of two waves v and w 
that propagate to the right and left, respectively, with speeds u + y/1 + 77 and u — •yTT^. 

Given a quasiuniform partition of [0,1] as in section 2, in addition to the spaces defined there, let for 
integer fc > 0 = {/ S 1] : /(I) = 0}, = {/ g 1),/(1) = 0}, and, for integer t- > 2, 


Sh — {(j) ^ C 


r-2 


\[Xj,Xj^x\ 


r_i, 1 < j < N}. Note that the analogs of the approximation and inverse 


properties (12.11) . (12.2L (12.51) . (12.6p hold for Sh as well, and that (12.3L (12.41) are also valid for the projection 

o 

V onto Sh, mutatis mutandis. 

The (standard) Galerkin semidiscretization of (ISW2al) is then defined as follows: Seek Vh ■ [0,T] —>■ Sh, 
Wh ■ [0, T] Sh, such that for t S [0, T] 


with 


{Vht,4>) + {{Uo + So)Vhx, 4>) + ^{VhVhx, 4>) + ^{whVhx, (l>) = 0, \/4>€ Sh, 

Q -1 O 

{wht,x) + ((wo - 5o)whx, x) + 2(whWhx, x) + 2^VhWhx, x) =0, Vx G Sh, 
VhiO) = P{v°), Whi0)=V{w°). 


(3.6) 

(3.7) 

(3.8) 


The main result of this section is 


ProDosition 3.1. Let (v.w) be the solution of llSW2a\) and assume that the huvotheses m and HY^) hold, 
that r > 3, and that h is sufficiently small. Then the semidiscrete ivp M-M has a unique solution 
{vh,Wh) for 0 <t <T that satisfies 


max ( X - Vh\ 
0<t<T 


+ IItx - TXTtll) < Ch 


r— 1 


If ( 77 , u) is the solution of (SWf^) and we define 


then 


Vh = [h{vh - Wh) + <5o]^ - 1, Uh = Vh + Wh + Uo, 


^max^dh- 77;i|l + ||7i-7ift||) < Ch^ 


(3.9) 


(3.10) 


Proof. Let p = v — Pv, 9 = Pv — Vh, a = w — Vw, ^ = Vw — Wh. After choosing bases for Sh and Sh we see 
that the ode ivp (I3.6l) - (l3.8p has a unique solution locally in time. From (ISW2al) and (13.61) . (13.71) we obtain, 
as long as the solution exists, 

Q 1 O 

{ 9 t , 4>) + ((mo + ^o){(ix + Px), 4>) + ^{vVx - VhVhx, 4>) + ffiiwvx - WhVhx, 4>) =0, yfi £ Sh, 

(6.x) + ((wo - So)i(Tx + f,x),x) + ^{wwx - WhWhx,x) + ^{uwx - VhWhx,x) = 0 . Vx G Sh. 
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(3.11) 

(3.12) 




























Now, since 


it follows that 


where 


VVx - VhVhx = {vp)x + {vO)x - {p0)x - PPx - 69x, 

WVx - WhVhx = w{px + Ox) + Vx{(T + 0 - {px + Ox){cr + ^), 
wwx - WhWhx = iwa)x + iw^)x - {crC)x - crax - ^^x, 
vWx - VhWhx = v{ax + ^x) + Wxip + 0) - {ax + ^x){p + 0), 

VVx - VhVhx = {v0)x - 00X + Rll, WVx - WhVhx = -0x^ + Rl2, 
WWx - WhWhx = {wCjx - ^£,x + i?21, VWx “ VhWhx = -^x0 + i?22, 


^11 = {vp)x - {p0)x - PPx, Ri 2 = WPx + W0X + Vxa + Vx^ - Pxa - Pxi - 0xa, 

R 21 = {wa)x - ia^)x - aax, R 22 = vax + v^x + WxP + Wx0 - axP - ax0 - ^xP- 

Putting now (f) = 0 ir^ (13.111) and x = ^ in (13.121) we obtain 

+((wo + ^o)0.,0) + U^v0)x,0) - l{09x,0) 

= - ((wo + (5o)px, 0) - |(i?ll, 0) + \{0xi, 0) - \{Ri2,0), 

hiur+iino - So)^x,o + §(«)., 0 - i(ee.,e) 

= - ((wo - (5o)o-a:,^) - |(i?21,C) + 5(?a:^iC) “ ^{R 22 , 0 ), 
Integration by parts yields (we suppress the t-dependence) 

{{uo + So)9x,0) = ^^^^0^1), i{v9)x,0) = livx0,9) + lv{l)0^1), 

{90X, 9) = ^03(1), ((wo - 5o)ix,0 = 

{{Ox,£,) = \{wxi,i) - \w{Q)e{^), {iix,i) = -\e{^)- 

Hence, (13.171) becomes 

\i\\0r + \{u, + 5, + lv{l)-9{l))0^{l) 

= - ((wo + So)px,0) - \{Vx0,0) + \ {0x^,0) - \ {Rll,0) - \ {Ri2,0). 
By (|Y2p and (13.5p we see that mo + 5o + |w(1) > cq >0. Therefore the above equation gives 

+ Kco - 0 ( 1 )) 0 ^( 1 ) < -((wo + 6 o)px, 9 ) 

- l{Vx9,0) + ^0x^,0) - U^11,0) - i(i?12,0). 

Similarly, from (j3.18p we obtain 

+ ^(-(«o - ^0 + fw;(0)) + m)ei0) 

= -{{uo - 6o)ax,0 + ^{^x9,0 - l{wx^,0 - §(-^ 21,0 - |(i?22,0- 
Again, by (IY2I) and (13.51) we get uq — ^o + |w(0) < —co < 0. We conclude that 

+ m)ei0) < -((wo - 

+ - 1(^21,0 - ^(^22,0- 
Finally, adding (I3.19P and p.20l) we get, as long as the solution of (I3.6I) - (I3.8I) exists, that 

kim\^ + iicf) + Kco - 0(i))0^(i) + i(co+ m)eio) 

< -iiuo + So)px,9) - ((wo - 6o)ax,0 - U'^x0,0) - j{wxi,0 
+ - 1 (^ 11 ,^) - 1 (^ 12 , 0 ) - |(i? 2 l ,6 - 5 (^ 22 , 0 - 


(3.13) 

(3.14) 

(3.15) 

(3.16) 


(3.17) 

(3.18) 


(3.19) 


(3.20) 


(3.21) 
















In view of (13.8p . by continuity we conclude that there exists a maximal temporal instance th > 0 such that 
Vh, Wh exist for t < th and 


||6'(i)||i.oo + Iie(i)lli.oo < Co, te [0,th]. 

(3.22) 

Suppose that th < !'. I’or t f \0.th \ we have bv p.22p 


i(co-0(l))0^(l) + Kco+C(O))f(O)>O, 

(3.23) 

and 


^i(0.e,0)i + ^i(e.0,e)i<|ii0|iiieii. 

(3.24) 

We obviously have 


I(^x0,0)i + i(«;.e,e)i<c^(ii0f+ iief). 

(3.25) 

Using now the approximation and inverse properties f|2.ip-(|2.6D for Sh (and also for Sh) 

we estimate the rest 


of the terms in the right-hand side of (13.211) as follows. We first clearly have 

|(K+<5oK,0)| + |((uo-<5oK,e)l +lieil)- (3.26) 


Integrating by parts we see by (13.151) that 

{Rii,0) = {{vp)^,e) - {{p9)^,0) - ippx,0) 

= - {vp, 0x) - p(l)d^(l) (p6», 0^) - ippx,0)- 

Therefore 

\{RlU0)\ < qiplloollelloo + C\\p\U\0x\\ + ||p||oo||0|lL + l|p||oo||0||||0.|| + ||p||oo||p.|l||0|| 

<Ch^\\0\\oo + Ch'^\\0x\\+Ch^\\0\\l + Ch^\\0\\\\0x\\+Ch^^-^\\0\\ (3.27) 

<ch^-\\\0\\ + \\0r). 

Integration by parts and (13.151) yield for the i?i 2 term that 

{Ri 2,0) = (wpx,0) - ^(wx0, 0) + {vxcr, 0) + 0) - (pxcr, 9) - {pxS,, 0) - {0xcr, 9). 

Hence, similarly as above 

\{Ri 2, 0)\ < Ch^-^\\0\\ + cpr + Ch^m + CIICII ||0|| 

+ Ch^'-^\\0\\ + C'/i^-'||?||oo||0|| + Ch^\\0\U\\0x\\ (3.28) 

<ch'^-^\\e\\ + c\\ 0 r + c\moi 

Again, using integration by parts and (13.151) for the i? 2 i term, we obtain 

(.^21,0 = -{wa,ix) - w( 0 )cr( 0 )^( 0 ) -h + cr( 0 )^^( 0 ) - (aaxA)- 

Therefore 

l(i?2i,e)l < qklllic.ll + qkiioolieiU + lklloo|iail?.ll + IklUlieilL + Iklloolk.imcil 

< ch^UxW + ch^UWoo + chnicilllC.II + ch^UWl. + Ch^^-^UW (3.29) 

<ch^-Hm\ + mn 

Finally, by (13.151) and integration by parts we have for the R 22 term 

{R22,0 = {vcrx,0 - ^{Vx^,0 + {WxP,0 + {Wx0,O - (.f^xP,0 - {<^x9,0 “ {pix,^- 

H 6 I 1 C 6 

\{R 22 , 0 \ < CIla.Illieil + Clief + CIIpIIIICII + C|| 0 ||||C|| 

+ Ik.||||p||oo||ell + Ik.||||0||oo||ell + I|p||oo||e.y^ 

< ch^-^uw+ teller+ ch^m+cm neii+( 3 . 30 ) 
+ ch’'-i||0|U||eil + c/inie.llll?ll 
<C'/1’'-I||cil + C||cf+ C||0||||ei|. 

By (13.211) . taking into account p.23l) - p.30l) we see that 

kitiPf + lief) < Ch^-\\\0\\ + lieil) + C{\\0f + ll^f), t e [0,th]- 

11 




















An application of Gronwall’s Lemma and (13.8p yield 

\\om + \\m\\<ch^~\ t€[oM. (3.31) 

from which by inverse assumptions it follows that ||0||i,oo + ||^||i,oo < for t G [0,th]. Since it was 

assumed that r > 3 this contradicts the maximality of th and (13.311) holds for 0 < t < T. The estimate (13.91) 
follows. Since now ||w — u?i||oo < IIpIIoo + Halloo < and similarly ||?ii — w/i||oo < and since 

V - Vh = [So + i((^ -w) + {vh - w/i))][(u -w)- {Vh - W/j)], 

we conclude that [[■q-'qhW < C{\\v-Vh\\ + ||w-w?i||) < Ch''~^. Similarly ||u-Uft,|| < Hu-u/iH + ||ic-w/iH < 
Ch^~^, and the proof of Proposition 3.1 is now complete. □ 

4. Numerical implementation and experiments 

4.1. Supercritical case. We will implement the standard Galerkin method for the shallow water equations 
with characteristic boundary conditions in the supercritical case using the space of piecewise linear continuous 
functions on a uniform mesh in [0,1] in the usual manner that problems with nonhomogeneous Dirichlet 
boundary conditions are approximated in practice. For this purpose we let Xi = ih, 0 < i < TV, Nh = 1, 
define Sh = {4> ^ C'°[0) 1] ^ ‘/’L i S Pi, 0 < j < TV — 1}, and let Sh consist of the functions in Sh that 
vanish at a: = 0. We seek r]h, uu '■ [O.T] —!> St, the semidiscrete approximation of the solution of (ISWII) . 
satisfying for 0 < t < T r]h{0,t) = rjo, Uh{0,t) = uq, and the system of ode’s 

o 

{r]ht,4‘) + {uhx,(/)) + {{r]hUh)x,4>) = y4'&Sh, , ^ 

(4.1) 

{Uht, 4>) + ivhx, (p) + {uhUhx, 4>) = 0> V()) G Sh, 

with initial values ?7;i(0) = Prf, Uh{0) = PvP, where P is the projection onto Sh- We discretize this ode 
ivp in time by the ‘classical’ explicit 4*^-order accurate Runge-Kutta scheme, written in the case of the ode 
y' = f{t, y), 0 <t <T, in the form 

y"4 = 2/” + |/(t” + |,y”), 

y”’" = 2/” + fc/(i” + fc,j/”’'), (4.2) 

= y- + fc(i/(t",J/”) + i/(r + |,J/"4) 

where k is the time step, t" = nk, n = 0,1,..., M — 1, Mk = T, and y" approximates y{t'^)- Theoretical 
and numerical evidence from linear stability theory and previous work by the authors, [Xm] . |Ad 2] , on 
similar nonlinear systems, suggests that the resulting fully discrete scheme is stable under a Gourant-number 
restriction of the form k/h < rg. 

In our first numerical experiment we check the spatial rate of convergence of this fully discrete scheme. 
We consider (ISWII) with lyo = 1 and uq = 3 and add right-hand sides to the pde’s so that the exact solution 
of the ibvp is r]{x,t) = xe~^* + ijo, u{x,t) = (1 — a; — cos(7ra;))e^* -I- ug- With h = 1/N, k = h/10 (so that 
the temporal error is negligible), we obtain the errors and associated rates of convergence of (essentially) 
the semidiscrete problem at T = 1 shown in Table 14.11 The experimental rates of convergence for both 


N 

V 

order 

u 

order 

40 

1.243098(-3) 

- 

5.623510(-3) 

- 

80 

3.110525(-4) 

1.99871 

1.405648(-3) 

2.00024 

160 

7.778520(-5) 

1.99959 

3.513979(-4) 

2.00006 

320 

1.944737(-5) 

1.99992 

8.784876(-5) 

2.00001 

480 

8.643341(-6) 

1.99998 

3.904381(-5) 

2.00001 

520 

7.364768(-6) 

1.99996 

3.326806(-5) 

2.00001 


Table 4.1. errors and spatial orders of convergence, supercritical case. 
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components of the solution are clearly equal to 2, i.e. superaccurate, as the expected rates for a general 
quasiuniform mesh would be equal to 1. A numerical study of stability for this example indicates that 
the errors remain of the same order of magnitude at T = 1 up to about k/h = 0.13. For larger values of 
k/h blow-up eventually occurs. It should be noted that for the same test problem the alternative standard 
Galerkin formulation (I2.7I) - (I2.9I) (analyzed in section 2) coupled with the same Runge-Kutta scheme gives 

errors that coincide with those of Table HIT] to at least 5 significant digits; the stability condition was also 
the same. (Recall that the result of Proposition 2.1 strictly holds for r > 3 due to the technical requirement 
in the proof for controlling the norm of the error. The numerical results suggest that the scheme 

converges for r = 2 as well and that the superaccurate order of convergence for a uniform mesh for r = 2, 
proved in [AD2) for an ibvp for the shallow water equations with homogeneous Dirichlet boundary conditions 
on u, persists in the case of the ibvp (ISWII) too.) 

Since the temporal error is much smaller than the spatial one, the experimental estimation of the temporal 
order of convergence may be done in the following way, used in |BDMK] . Let be the fully discrete 


k/h 

E*{T) 

order 

E{T) 

1/35 

2.6618459890(-8) 

- 

1.9910684230(-4) 

1/40 

1.6020860073(-8) 

3.8022 

1.9910680992(-4) 

1/45 

1.0112973048(-8) 

3.9061 

1.9910679532(-4) 

1/50 

6.6717108025(-9) 

3.9478 

1.9910678792(-4) 

1/55 

4.5726218272(-9) 

3.9638 

1.9910678370(-4) 

1/60 

3.2362144361(-9) 

3.9728 

1.9910678102(-4) 

1/64 

2.5020819256(-9) 

3.9865 

1.9910677950(-4) 

1/64.5 

2.4254282105(-9) 

3.9983 

1.9910677934(-4) 

1/65 

2.3516195603(-9) 

4.0020 

1.9910677918(-4) 


Table 4.2. Temporal order of convergence, supercritical case, scheme (|11])-(|12]), h = 
1/100, T = 1, kref = h/UO. 


approximation of For a fixed value of h we make a reference computation with a very small value k = 

kref - The approximate solution iL™ = H™{h, kref), where mkref = T, differs from the exact solution ? 7 (-, T) 
by an amount which is practically the error of the spatial discretization. For the same value of h we define a 
modified error for small values of k, that are nevertheless considerably larger than kref, by the formula 
E*{T) = \\H^{h,k) — H]^{h,kref)\\, where nk = T. Since taking the difference iL))(/i, A:) — HJ^{h,kref) 
essentially cancels the spatial error of HJ^{h,k), we expect that E*(T) will decrease at the temporal order 
of convergence of the scheme as k decreases. This is illustrated in the case of the test problem under 
consideration and the fully discrete scheme (|4. Ill - (14.21) in Table IT^ where h = 1/100, T = 1, kref = h/\20, 
and E{T) denotes the LF' error \\H'^{h,k) — ?7(t")||. For this range of fc’s the expected temporal order of 
convergence, equal to 4, clearly emerges. The analogous experiment with the Galerkin method (|2.7I) - (I2.9I) 
discretized in time with the same Runge-Kutta scheme yields fourth-order temporal convergence in as 
well. 

In the next numerical experiment we integrate (ISWlj) with the fully discrete scheme (j4.1l) - (l4.2|) . taking 
h = 1/N, N = 2000, k = h/10, rjo = 1, uq = 3, and initial conditions ry°(a;) = 0.05 exp(—400(a; — 0.5)^) -I- ryo, 
u°{x) = 0.1 exp(—400(x — 0.5)^) -I- mq) 0 < x < 1. (Small-amplitude initial conditions were taken to ensure 
that no discontinuities in the derivatives of the solution develop before the wave profiles exit the spatial 
interval of integration.) The evolution of the numerical solution is depicted in Figure Id.ll aVffL (The 
approximate ry-profiles are on the left and those of u on the right.) The initial Gaussian perturbations of the 
uniform state ?yo = 1, mq = 3 evolve into two unequal pulses for both ry and u that travel to the right and 
exit the computational domain by about t = 0.4 without leaving any visible residue or backwards-travelling 
oscillations as is confirmed by Fig. I4.11 gl that shows the time history of the quantities ma.Xx\ri{x,t) — 7yo| 
and maxa;|u(x, t) — uo|. (Here ?y, u denote the approximate solution.) Due to the presence of the spatial and 
temporal discretizations, the numerical boundary conditions are not expected to be exactly transparent. 
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u at t = 0.05 
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(d) T] and u aX t = 0.15 




(e) T] and u at t = 0.3 




(f) 77 and u at t = 0.4 
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(g) jn&yix\ri(x,t) — ? 7 o| and iiiexx\u{x,t) — uo| vs. time 




(h) Magnification of (g) 


1.67 

1.66 

1.65 

1.64 

1.63 

1.62 

1.61 

1.6 

1.59 

1.58 



(i) maxa;(M— y/1 + rj) vs. time 

Figure 4.1. Evolution of Gaussian initial profiles, supercritical case. 


However, they are highly absorbing; Figure HTlT hl reveals that the residue after the waves exit is of 0(10“®). 
The positivity of Tn&yix{u — + 77 ) for all t checked in Fig id.lf il confirms that the numerical solution has 

remained supercritical throughout the evolution. 

In order to study numerically the stability of the fully discrete scheme for this test problem, 

as there are no exact solutions, we took as a measure of error the residual quantity max^, I 77 — ryol- This is 
plotted in Figure|4Tjg) and (h), stabilizes after the waves exit the computational domain, and has the value 
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9.76i? — 07 at t = 0.45. We then increased k and observed that this residual was conserved up to about 
k/h = 0.3695 and started increasing afterwards. Since the maximum wave speed c is the speed of the higher 
rightward-travelling pulse, which is equal to u -I- -|- rj ~ 4.5 for the duration of this experiment, we obtain 

a Courant number restriction of about ck/h < 1.67. (Linear stability theory for this method applied to the 
model problem rjt -I- cijx = 0 with periodic boundary conditions and c constant would give a Courant number 
restriction ck/h < -^8/3 ~ 1.633 which is not far from the experimental result for this small-amplitude 
nonlinear propagation problem.) It should be noted that the scheme (I2.7I) - (I2.9I) discretized in time with the 
Runge-Kutta method (14.21) yields practically the same numerical results for this test problem. 

4.2. Subcritical case. We implement the standard Galerkin method for the SW with characteristic bound¬ 
ary conditions in the subcritical case using again piecewise linear continuous functions on a uniform mesh in 
[ 0 , 1 ]. In addition to the notation introduced in the previous subsection for the mesh on [ 0 , 1 ] and the space 
Sh, we let Sh,o consist of the functions in Sh that vanish at cc = 0 and x = 1. We first consider the ‘direct’ 
standard Galerkin semidiscretization of (ISW2|) . i.e. without reducing first the ibvp into the diagonal form 
(ISW2a)) . We seek accordingly u/j : [0, T] —>• the semidiscrete approximation of the solution of (ISW2I) . 

satisfying for 0 < t < T 

{r]htA) + {uhxA) + {{'nhUh)xA) y4>&Sh, (4.3) 

{Uht, x) + ivhx, x) + (uhUhx, x) =0, Vy e Sh,o, (4.4) 

where Uh G Sh,o, and Uh{xi,t) = Uh{xi,t), I <i < N - I, 

Uh{xQ,t) = -2sjl + rih{xQ,t) +Uo + 2y^l + r]o, (4.5) 

UhixN,t) = 2^/Tl^ri/J/x//^ + uo-2^/TT^. (4.6) 

At t = 0 we compute Uh{0) and r]h{0) as the projections of the initial data 77 ° , onto Sh- Although the 
semidiscrete ivp (14.311 - (14.61) has nonlinear boundary conditions, it is easily discretized in time by an explicit 
scheme such as the 4*^-order RK method ()4.2L by first advancing from t" to the approximations of rjh 
and of the ‘interior’ Uh using the temporal discretizations of (j4.3l) and (14.4L and then updating the Uh{xQ, t) 
and Uh{xN,t) values at t = by (14.5p and (14.6p . 

The convergence of this scheme was not analyzed in section 3. The following experiment suggests that 
in the case of uniform mesh its errors are of 0{h^). We consider (|SW2p with rjo = 1, uq = 1 and its 


N 

V 

order 

u 

order 

40 

4.847892(-3) 

- 

2.932354(-3) 

- 

80 

1.207564(-3) 

2.00526 

7.414336(-4) 

1.98367 

160 

3.017313(-4) 

2.00076 

1.860285(-4) 

1.99479 

320 

7.544641(-5) 

1.99974 

4.657627(-5) 

1.99786 

480 

3.353298(-5) 

1.99991 

2.071174(-5) 

1.99867 

520 

2.857355(-5) 

1.99953 

1.764866(-5) 

1.99944 


Table 4.3. errors and spatial orders of convergence, subcritical case, semidiscretization 
(IdAll-dTbll. 


semidiscretization (IQli-dreli with piecewise linear continuous functions. We add appropriate right-hand sides 
to the pde’s in (ISW2I1 so that the exact solution of the ibvp is r]{x, t) = {x + l)e~^*, u{x, t) = (2x-|-cos(7ra;) — 
l)e*+xA{t) + {l-x)B{t),wh.eTe A{t) = 2y/T^^//(/r/t)+uo-2^/Tl^, B{t) = -2yT+"^^0)^-|-uo-l-2vTT^. 
We consider uniform spatial and temporal meshes with h = 1/N, k = h/10, and discretize the semidiscrete 
problem in time using again the 4*^-order ‘classical’ RK scheme. (We checked that the temporal error 
is negligible for the range of N’s that we tried.) The resulting errors and rates of convergence of 
(essentially) the semidiscrete problem at T = 1 are shown in Table H751 The rates are practically equal to 2, 
i.e. superaccurate, as in the supercritical case. An analogous temporal-order calculation to that appearing in 
Table [4?^ was not so robust and gave rates between 3.7 and 3.9; thus some sort of temporal order reduction 
cannot be ruled out for this scheme. The experimental Courant number restriction was k/h < 0.53. 
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In the following numerical experiment we integrate (ISW2|) with the same fully discrete scheme, taking 
h = 1/iV, N = 2000, k = /i/lO, rjo = uq = 1, and initial conditions r]^{x) — 0.1 exp(—400(x — 0.5)^) + 770 , 
u°(x) = 0.05 exp(—400(x — 0.5)^) +uo. The evolution of the numerical solution is shown in Figure IT^ lal-fil. 




(a) T] and u at t = 0.0 
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(b) Tj and u at t = 0.05 
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(h) rj and u a.t t = 1.25 
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(j) rj and u at t = 1.5 




(k) Tiiayix\ri{x,t) — ryo] and Tnayix\u{x,t) — mo| vs . time 














(1) manification of (k) 



(m) maxa;(u — ^1 + rj) vs. time 

Figure 4.2. Evolution of Gaussian initial profiles, subcritical case, semidiscretization (1431)-(lEGl). 

The initial Gaussian perturbations of the steady state ryo = 1, ito = 3 evolve into two unequal pulses for 
both components of the solution, which travel to the right and left and exit the computational domain at 
about t = 0.25 and t = 1.35, respectively, without leaving any visible residue as confirmed by the temporal 
history of the maximum deviations of the approximations of ry and u from ryo and uq shown in Figures l4.2f kl 
and (1). The latter graph shows that the maximum residue after exit is of 0(10“®), confirming the high 
degree of absorption of the discrete characteristic boundary conditions. The quantity maxa;(it — i/l -|- rj) 
remains negative, i.e. the numerical solution is subcritical, throughout the evolution, cf. Fig. \4:.2{ m). We 
investigated the stability of this fully discrete scheme by using again as a measure of error the quantity 
maxa; |ry — 7yo| which stabilizes at t = 1.55 to the value 1.21F — 05. The maximum wave speed c for this 
problem is about 2.5, and the observed Gourant number restriction was ck/h < 1.67 in conformity with the 
analogous value in the supercritical case. 

In section 3 we analyzed a different Galerkin semidiscretization of the ibvp under consideration, namely 
that given by (IT61)-(I331). a semidiscrete approximation of the diagonal form of the ibvp, i.e. of (ISW2al) . As 
the following experiment suggests, this semidiscretization is also 0{h?) accurate in if we use piecewise 
linear continuous functions on a uniform mesh. We consider the inhomogeneous version of (ISW2al) with 
unknowns 2v and 2w instead of v and w and take as exact solution the functions v = u — uo + 2{^/l + rj — So), 
w = u — uq — 2{y/n^ — (5o), where So = ^1 -|- ryo and r]{x,t) = {x -\- u{x,t) = (2x + cos(7ra:) — 

l)e* -I- xA(t) -b (1 — x)B{t), with A{t) = 2-^/1 + ry(l,t) -b mq — 2(5o, B{t) = —2\/T+ij{0^ -b uo + 2(5o. For 
r]o = Uq = 1 we approximate this ibvp by the nonhomogeneous analog of the semidiscretization (I3.6I) - (I3.8I) 
(with unknowns 2vh and 2wh) using piecewise linear continuous elements with h = 1/N, and the ‘classical’ 
4*^-order RK scheme for time stepping with k = h/10. The resulting errors and rates of convergence 
at T = 1 are given in Table 14.41 The rates are practically equal to 2 as in the previous cases, due to the 
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N 

1 

order 

u 

order 

40 

2.470369(-3) 

- 

9.918820(-4) 

- 

80 

6.172661(-4) 

2.00076 

2.472869(-4) 

2.00398 

160 

1.543038(-4) 

2.00012 

6.179903(-5) 

2.00053 

320 

3.857665(-5) 

1.99997 

1.545737(-5) 

1.99929 

480 

1.714531(-5) 

1.99998 

6.870865(-6) 

1.99967 

520 

1.460903(-5) 

1.99999 

5.854663(-6) 

1.99958 


Table 4 . 4 . errors and spatial orders of convergence, subcritical case, semidiscretization 

(EHi-dall). 


uniform mesh. The temporal order of convergence for this scheme was found to be practically equal to 4; 
the associated temporal-order calculation with h = 1/50, kref = h/200 was quite robust. 

We now compare the solutions of the two semidiscretizations (ig-(ITO and dSlli-dH by means of 
a numerical experiment. We consider the ibvp (ISW2I) with rj^ = uq = 1 and initial values ri(x, 0) = 
0.1 exp(—400(a; — 0.5)^) u{x, 0) = 0.05 exp(—400(a:: — 0.5)^) (Recall that the temporal evolution of 
the numerical solution of this ibvp generated by (|4.3I) - (I4.6I1 with N = 2000, h = 1/A^, k = h/lO is shown in 
Figure 1121) We will compare the numerical solutions of this ibvp with both discretizations, using the same 
numerical parameters. Let [rjh, Uh) denote the fully discrete approximations at time t produced by (I4.3I) - (I4.6I) 
as underlying semidiscretization. In addition, let (rihD^Uho) be functions in Sh with point values computed 
by the formulas (13.1011 . where Vh, Wh are now the fully discrete approximations at t produced when we use 


method (|3.6ll-(|3.8ll. 

Define e(t) = 

maxo<i<Ar 

\r]h{xi,t) - 

dhoixi^t)] and 

e(t) = maXo<i<Ar M?i(: 

t 

0.0 

0.05 

0.1 

0.15 

0.2 

0.225 

0.25 

e{t) 

2.0(-8) 

2.0(-8) 

2.0(-8) 

2.0(-8) 

4.0(-8) 

3.0(-8) 

3.0(-8) 

e{t) 

< 10“® 

< 10“® 

< 10“® 

< 10“® 

2.0(-8) 

2.0(-8) 

2.0(-8) 


t 

0.35 

0.75 

1.05 

1.15 

1.25 

1.35 

1.5 

e(t) 

1.232(-5) 

1.301(-5) 

1.237(-5) 

1.3(-5) 

1.224(-5) 

1.23(-5) 

1.217(-5) 

e{t) 

8.71(-6) 

9.2(-6) 

8.75(-6) 

9.54(-6) 

8.65(-6) 

8.75(-6) 

8.6(-6) 


Table 4.5. Comparison of the discretizations (14.311 - (14.611 and (I3.6I1 - (I3.8I1 . subcritical case, 
experiment of Fig. 14.21 


UhD{xi,t)\. The quantities e and e for various values of t = f" G [0,1] are given in Table 021 We observe 
that up to about t = 0.25 the errors are of 0(10“®) and subsequently increase to 0(10“^). This is due to 
the fact that about t = 0.25 the larger, rightwards-travelling wave completes its exit from the computational 
domain (see Figure 02Ke), and we expect a small residue to be radiated into [0,1] because the numerical 
characteristic boundary conditions of (I4.3I1 - (I4.6I1 are not exactly transparent. This residue has a magnitude 
of 0(10“®) as evidenced by Fig. I4.2l lb It is worthwhile to note that the ‘diagonal’ approximation (I3.6II - (I3.8I1 
leaves a practically negligible residue. This is suggested by the evidence in Figure 021 In this figure, the two 
graphs on the left depict what is left in the computational domain of the ry-and M-components of the numerical 
solution generated by (14.31) - (14.61) at t = 1.5, after the waves have exited the domain, while those on the right 
are the analogous ry-and M-profiles generated by (13.611 - p. 81) . (Thus the graphs on the left are magnifications 
of the graphs of Fig. mi)-) We observe that the residues of the usual Galerkin semidiscretizations are of 
0(10“®) while those of the ‘diagonal’ scheme are much smaller. It is clear, at least for this example, that 
the ‘diagonal’ Galerkin method has practically transparent boundary conditions. 
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(a) rjh and -quo, i = 1.5 



Figure 4.3. Magnifications of the profiles of the numerical solutions generated by the 
methods (14.31) - (14.61) (left) and (I3.6I) - (I3.8I) (right) at t = 1.5. Subcritical case, evolution as in 
Fig. 14.21 


4.3. Linearized vs. nonlinear characteristic boundary conditions in the subcritical case. If the 

elevation of the free surface ry is a small perturbation of the steady state 770 one may derive linearized 
approximations to the characteristic boundary conditions: Assume that the wave height is given by 1 -|- ?y = 


l + ?7o + i?, where |? 7 | < 1 - 1 - 770 . Then ^l + q = -v/1 + ?7o + : 


n-no 


2\/1-}-77o 

from the characteristic boundary condition, for example at a; = 1, it follows that 


0{rp), and, in the subcritical case (ISW2I) . 


u{l,t) - 


VI + Vo 


= Uo 


Vo 


VI + Vo 




from which neglecting the 0(rj^) term we obtain the linearized b.c. 


and similarly at a: = 0 the b.c. 


u{l,t) - 




^(ly) 

VI + Vo 

v{o,t) 


= Uo - 


= Uo + 


Vo 


V^ + Vo’ 


Vo 


(4.7) 


(4.8) 


Vo V^ + Vo' 

The linearized boundary conditions (I4T1) and (j4.8|) have often been used in the computational fluid 
dynamics literature, as was mentioned in the Introduction. It is not hard to see that (1471) and (ra are 
transparent for the linearized shallow water equations obtained by linearizing the system about the 
steady state {uo,Vo)- This may be seen by diagonalizing the linearized system and explicitly computing its 
solutions propagating along the characteristics. In [SLTT] it is proved that the ‘energy’ integral /giu^ + 
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Y^;^rf)dx of the solution (? 7 ,u) of the linearized system decreases with time in the presence of a class of 
boundary conditions that includes (I4.7p - (j4.8ll . a fact that implies the well-posedness of the ibvp for the 
linearized system supplemented by Kn and (ira . 

However, the linearized characteristic boundary conditions are not transparent for the nonlinear system 
(HU). So, there arises a need to study their absorption properties. In [NMF] and [SLTT] the linearized 
conditions were compared with the nonlinear ones, appropriately discretized in the finite difference schemes 
used in these two references. The comparison was effected by means of numerical experiments, many of which 
in more complicated instances of hydraulic and geophysical interest, including single- and multi- layered flows 
in the presence of variable bottom topography, cross velocity terms, and Coriolis forces. Moreover, the finite 
difference schemes used in |NMF) and |SLTT| allow in general the simulation of flows that develop steep 
fronts and discontinuities. It was confirmed that the linearized conditions in many examples give rise to 
spurious oscillations that are reflected backwards into the computational domain. 

In the sequel we will compare the two sets of boundary conditions in the case of smooth subcritical flows 
discretized by the ‘direct’ fully discrete Galerkin method with piecewise linear continuous functions in space 
coupled with ‘classical’ 4‘^-order RK time stepping. The spatial discretization in the nonlinear b.c. case is 
given by (ig-dTBl): in the linearized case gU-gU) are replaced by their obvious linearized analogs resulting 
from gU-gU)- We first check the errors and associated orders of convergence for the scheme with the 
linearized b.c.’s. Adding appropriate right-hand sides to the pde’s in (ISW2I1 so that the exact solution of the 
ibvp is T](x, t) = {x + I)e“‘^*, u{x, t) = {2x + cos(7ra:) — l)e‘ -I- xa{t) -I- (1 — x)b{t), where a{t) = uq -I- ’ 

b{t) = uo + Ji+lg ’ uo = rjo = I and compute with h = 1/N and k = h/20. The resulting errors 

of (essentially) the spatial discretization are shown, accompanied by the associated orders of convergence in 
Table lT6l Due to the spatial uniform mesh the rates are again practically equal to 2. 


N 

1 

order 

u 

order 

40 

4.835002(-3) 


2.930984(-3) 


80 

1.204245(-3) 

2.00539 

7.408500(-4) 

1.98413 

120 

5.350289(-4) 

2.00088 

3.299726(-4) 

1.99472 

160 

3.008683(-4) 

2.00099 

1.858783(-4) 

1.99497 

200 

1.925597(-4) 

1.99991 

1.190431(-4) 

1.99695 

240 

1.337304(-4) 

1.99966 

8.269369(-5) 

1.99835 

280 

9.825114(-5) 

1.99998 

6.077487(-5) 

1.99783 

320 

7.523153(-5) 

1.99920 

4.653473(-5) 

1.99936 

360 

5.944094(-5) 

2.00018 

3.677599(-5) 

1.99820 

400 

4.814899(-5) 

1.99964 

2.979144(-5) 

1.99908 

440 

3.979234(-5) 

2.00006 

2.462384(-5) 

1.99880 

480 

3.343735(-5) 

1.99975 

2.069270(-5) 

1.99898 

520 

2.849223(-5) 

1.99946 

1.763197(-5) 

1.99978 


Table 4.6. errors and spatial orders of convergence, subcritical case, semidiscretization 
(14.31) . (14.41) with the linearized b.c. (14.71) . (14.81) . 


We now consider a numerical experiment discussed in Section 4.2, wherein the fully discrete numerical 
scheme based on the semidiscretization (14.31) - (14.61) gives the evolution depicted in Figure HU We repeat 
the experiment with the linearized boundary conditions gu, gu instead of (I4.5L(I4.6I) . In Figure Hg] 
the graphs of rjun, uun (solid line), i.e. the numerical solution corresponding to the linearized b.c.’s, are 
superimposed on those of the numerical solution unuu (dotted line) computed with the nonlinear 

b.c.’s. At t = 0.25, when the rightwards-travelling pulse has almost completed its exit from the computational 
domain at x = 1 (recall Fig. l4.2l gH. we observe that the linearized b.c.’s give rise to a spurious pulse of 
amplitude about 4 x 10“^ that is reflected inwards and travels (as Fig. Id.df bl confirms) to the left with a 
speed of —uq -\- + rjQ = 0.41 as predicted by the analysis of the linearized system of which this small- 
amplitude pulse is an approximate solution. A similar (but much smaller) rightwards-travelling pulse is 
created at x = 0 when the smaller wave in Fig. 14.21 exits the computational interval. (These pulses will 
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(a) T]iin. VS. rjNiin, and uun. vs. umuu. at i = 0.25 




(b) rjiin. vs. r]Miin. and uun. vs. mathk. at t = 0.35 

Figure 4.4. Reflected spurious pulses due to the linearized b.c.’s (solid lines) superimposed 
on the solution of (14.31) - (14.61) (dotted lines); evolution of Fig. 14.21 


eventually exit the interval since the linearized b.c.’s are transparent for the linearized system.) This may 
also be seen in the similar numerical experiment shown in Figure H3] We consider (ISW2I) with r/g = uq = 0 




(a) ryiin. vs. ’qMUn. and uun. vs. unuu. at t = 0.6 
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0.001 




(e) rjiin. vs. ■qNlin. and UUn. vs. UNlin. at t = 1.6 

Figure 4.5. Reflected spurious pulses due to the linearized b.c.’s (solid lines) superimposed 
on the solution of (I4.3|) - (l4.6p (dotted lines); evolution resulting from rf = 0.1 exp(—400(a;— 
0.5)^), = 0, ryo = Mo = 0. 


and initial conditions rf{x) = 0.1 exp(—400(x —0.5)^), vP{x) = 0, that we solve by the fully discrete Galerkin 
method corresponding to (I4.3l) - (l4.4p with nonlinear and linearized b.c.’s using h = 1/N, k = h/10, N = 2000. 
Due to symmetry the initial Gaussian r]^(x) breaks up into two waves that travel with equal speeds in opposite 
directions and exit the computational domain at about t = 0.6. Figure lDSl show magnifications of the ensuing 
behavior of the numerical solution computed with the linearized b.c.’s (solid line) superimposed on the correct 
steady state (dotted line), the result of computing with the nonlinear b.c.’s. Two equal spurious pulses are 
reflected inwards at both ends and, according to linear theory, travel with unit speed, interact linearly, and 
exit cleanly by t = 1.6. 

We consider next the shallow water equations written in dimensional variables {h,u) where u is the 
horizontal velocity and h denotes now the height of the water column above the bottom; the latter is located 
at a depth H below the level of rest. The system, posed on a channel of length 2L, is written, in the notation 
of [NMF) . in the form 


ht + {hu)x = 0, 
ut + ghx + uux = 0 , 


- L<x<L, t>0, 


(4.9) 


with initial conditions 

h{x,0) = f{x), u(x,0)=v(x), —L<x<L, (4-10) 


where g is the acceleration of gravity. The system is supplemented by the nonlinear characteristic boundary 
conditions now written as 


u{-L, t) + 2^/gh{-L,t) = a'^, (4.11) 

u{L, t) - 2^/gh{L,t) = a^, (4.12) 

where = uq ± 2i/(y/io, and uq, hg are the constant values of u and h outside [—L, L], i.e. the ‘asymptotic’ 
state of the flow. The linearized boundary conditions (obtained by linearizing (j4.11l) - (j4.12l) assuming that 
h = H + h, where h is small) are 


u{-L, t) + ^ j^h{-L, t) = b%, 
u{L,t) - ^J^h{L,t) = b+, 


(4.13) 

(4.14) 


where 5^ = uq ± y^hg. We repeat the numerical experiment in section 4.1 of |NMFj . using now the fully 
discrete Galerkin method with piecewise linear continuous functions in space and the 4*^-order classical RK 
method in time, and taking L = 2m., H = 0.2m, g = 9.8m/sec^, ho = 0.2m, uo = 0, i.e. = ±2.8m/sec, 
6^ = ±1.4 m/sec, and initial values /(x) = 0.25 m and v{x) = 0, with discretization parameters Ax = 4/7V, 
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k = At — 1/(10A^), N = 8000. We first solve the problem with the nonlinear b.c.’s (14.111) . (14.121) . The 
evolution of the waveheight h is depicted in Fig. 4.6. Due to the difference of the initial profile h{x, 0) = f{x) 
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0.05 


-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 


t = 2.0 


t = 3.0 


Figure 4.6. Evolution of waveheight h{x, t); numerical solution of (I4.9I) - (I4.12I) with f{x) = 
0.25 m, v{x) = 0, Hq = 0.2 m, uq = 0. 
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Figure 4.7. Graphs of and vs. time. Evolution of Figure H751 nonlinear 

and linearized characteristic b.c.’s 


at a; = ±2 with the asymptotic state ho, fluid exits the domain from both boundary points and two wavefronts 
are created that travel in opposite directions, interact, and completely exit the computational domain at 
about t = 2.25 sec with no apparent reflections. After this time the asymptotic state hg = 0.2 m is achieved 
throughout the domain. The time history of the solution {h,u) at a; = —1 is given in Figure [4.71 and 
corresponds to the three stages of the evolution deduced from Figure 14.61 The graph of the numerical 
solution corresponding to the linearized b.c.’s (14.13^ - 114.141) is superimposed on the previous graph. It may 
be seen that solution with the linearized b.c.’s is just slightly less absorbing. These results are in agreement 
with those of [NMFj . 

In our final numerical experiment we solve again the SW on [—L, L] in their dimensional form (14.9L (14.101) 
with the nonlinear b.c.’s (I4.11L (14.121) and the linearized ones (14.131) . (14.141) . taking now L = Im, H = 0.2m, 
g = 9.8m/sec^, ho = 0.2m, uq = 0, i.e. = ±2.8m/sec, 6^ = ±1.4m/sec, and as initial values in (|4.10l) 
v{x) = 0 and f{x) given by a half-sine pulse centered at x = 0 and equal to 0.2-1-0.05 sin(7r(x-|-0.3)/0/6) m if 
|x| < 0.3m and to 0.2m for |x| > 0.3m. The mesh parameters are Ax = 2/A, N = 4000 and At = 1/(10A). 
Figure shows the ensuing evolution of the numerical solution {h is shown on the left and u on the right 
at each temporal instance.) The graphs corresponding to the nonlinear characteristic b.c.’s (dotted line) are 
superimposed on those computed with the linearized b.c.’s. Two pulses are produced (with no oscillations 
apparent at the points where the spatial derivative is discontinuous) that travel to opposite directions and 
exit the computational interval [—1,1] with no reflections in the case of the nonlinear b.c.’s, shortly after 
t = 0.8 sec. The linearized b.c.’s give spurious reflected pulses with waveheights of amplitude of 




(a) h and u at t = 0.2 
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(k) h and u a.t t = 2.4 

Figure 4.8. Evolution of h{x, t) (left) and u(x, t) (right); numerical solution of (I4.9I) - (I4.10I) 
on [—1,1] with nonlinear (dotted line) and linearized (solid line) outflow b.c.’s with v{x) = 0, 

/(x) = 0.2 + 0.05 sin(7r(x + 0.3)/0.6) if |x| < 0.3, /(x) = 0.2 if |x| > 0.3. 

that travel and interact approximately linearly until they exit the domain at about t = 2.4 sec. The spurious 
reflections may also be observed in the graphs of the temporal history of the waveheight at a gauge at 
X = —0.5 m computed with the nonlinear (dotted line) and the linearized (solid line) boundary conditions. 




(a) h at X = —0.5 vs. time (b) magnification of (a) 

Figure 4.9. Graphs of h{—0.5,t) for 0 < < < 2.5; evolution of Figure iTSl Nonlinear b.c.’s 

(dotted line) vs. linearized b.c.’s (solid line). 

The passage of the main pulse at the gauge location is completed at about t = 0.6 sec. Subsequently the 
rightwards-travelling reflected spurious small-amplitude pulse due to the linearized b.c.’s is observed to be 
passing the gauge location at around t = \ sec, while the leftwards-travelling companion spurious pulse passes 
at around t = 1.7 sec. 

This experiment is qualitatively analogous to that described in Section 4.2 of [NMFj . We simulated here a 
smaller-amplitude wave to avoid the emergence of discontinuities during the temporal interval of the evolution 
of Fig. 14.81 In [NMF] the analogous maximum waveheight was equal to about 0.3 as the finite-difference 
scheme used in that work had artihcial viscosity (private communication by Dr. A. McC. Hogg), that 
smoothed out the emerging discontinuity. (A more correct simulation of the numerical experiment in [NMF] 
requires that an initial half-sinusoidal impulse is imparted on the water column. This may be modelled by 
an appropriate forcing term of the form —F{x,t)/h on the right-hand side of the momentum equation, i.e. 
the second pde in (14.91) . where F{x, t) may be taken as a multiple of sin(7rt) if (x, t) G [—0.1, 0.1] x [0,1] and 
F = 0 otherwise, while the initial conditions should now be /(x) = 0.2 m and v{x) = 0. For small enough 
amplitude of the forcing term F we observed that a qualitatively similar evolution took place.) 
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